reg.mean=function(data,x,y,N,mux){
  data=as.name(data); x=x; y=y; n=length(x); N=N; mux=mux
  fpc=(1-n/N); fit=lm(y~x); ybar=mean(y)
  plot(x,y,main=data); abline(fit); xbar=mean(x)
  b=coefficients(fit)[[2]]; SSE=deviance(fit)
  dfe=df.residual(fit); MSE=SSE/dfe
  muhat.yL=ybar+b*(mux-xbar)
  vhat.muhat.yL=fpc*(MSE/n)
  bound=2*sqrt(vhat.muhat.yL)
  lower=muhat.yL-bound; upper=muhat.yL+bound
  corr=cor(x,y)
  res=rstudent(fit); pred=fitted(fit)
  hist(res); plot(pred,res,main='Residuals vs. Predicted'); abline(0,0)
  qqnorm(res); qqline(res)
  cat("","\n","Results from SRS: Data =",data,"\n",'Estimation method = Regression for mean y','\n',"N =",N,"n =",n,
      "FPC =",fpc,"\n",'Slope b =',b,'','mux =',mux,'\n',
      "MuyL =",muhat.yL,"\n",'Correlation =',corr,'\n','Mean square error =',MSE,'\n',"Vhat muhat.yL =",vhat.muhat.yL,
      "\n","Bound =",bound,"\n","Lower Bound =",lower,"Upper Bound =",upper,"\n","")
  results=list(mux=mux,b=b,data=data,n=n,N=N,fpc=fpc,x=x,y=y,muhat.yL=muhat.yL,vhat.muhat.yL=muhat.yL,bound=bound,
               MSE=MSE,dfe=dfe,lower=lower,upper=upper,corr=corr,res=res,pred=pred)  
}

# to use the function with its call:
# reg.mean(data,x,y,N,mux)

# data: name of dataset, in quotes
# x: numeric vector (data)
# y: numeric vector (data)
# N: population size
# mux: mean of X